ccc

input data

load Epi and mobility data

d <- readRDS(file='Rdata/input/data_2020-04-12.rds')
#
country <- as.character(d$Country)
# mean sd of ditrivution infection to death
mu_id = 18.8
si_std = .45*mu_id

load and clean mobility data.

##  [1] "Algeria"                  "China"                   
##  [3] "Czechia"                  "Dominican_Republic"      
##  [5] "Ecuador"                  "Iran"                    
##  [7] "Peru"                     "South_Korea"             
##  [9] "United_Kingdom"           "United_States_of_America"
## [1] "warning"
## [1] "Algeria"            "China"              "Dominican_Republic"
## [4] "Ecuador"            "Iran"               "Peru"              
## [7] "South_Korea"

useful data

# correct epi
D <- d$D_active_transmission
mD <- as.matrix(D[,-1])
I <- d$I_active_transmission

cumD <- data.frame(D$dates,apply(D[,-1],2,cumsum))
cumI <- data.frame(I$dates,apply(I[,-1],2,cumsum))

N_geo <- ncol(D)-1
N_days <- nrow(D)
country <- as.character(d$Country)

process mobility data

load data (dots), get weekly average, interpolate between (blue line)

process Epi data

load Epi data and serial interval

delay infection to death

## [1] 0.0292052

Useful formatting

see paper

# delay death infection matrix
H <- matrix(0,nrow = N_days, ncol = N_days)
for (i in 1:N_days){
  f <- max(c(1,i-delta_id$SItrunc))
  H[i,f:i] <- rev(delta_id$dist)[((delta_id$SItrunc+1)-(i-f)):(delta_id$SItrunc+1)]
  if (i>1) H[i,f:i] <-  H[i,f:i]/sum( H[i,f:i])
}

# delay SI death to death matrix
W <- matrix(0,nrow = N_days, ncol = N_days)
for (i in 1:N_days){
  f <- max(c(1,i-SI$SItrunc))
  W[i,f:i] <- rev(SI$dist)[((SI$SItrunc+1)-(i-f)):(SI$SItrunc+1)]
  # if (i>1) H[i,f:i] <-  H[i,f:i]/sum( H[i,f:i])
}

# overall infectivity matrix
Ot <- W%*%mD
# correct infectivity when >0 case observed but no infectivity
for(i in 1:N_geo){
  Ot[ which( (D[,i+1]>0) & (Ot[,i] ==0 ) ),i] <-NA
}



# matrix of mobility and delayed version
M <- as.matrix(mobility$mob_combined_smooth[,-1])
# M_D <- H %*% M

Baseline model: estimate Rt

get epiestim Rt estimates with 7 day time window.

full model

parameters

# epi
R0 <- 5
# for risk
b <- 0.5

theta0 <- c(rep(R0,N_geo),
            rep(b,N_geo))

prior_theta <- matrix(c(rep(c(0,10),N_geo),
                        rep(c(0,1e3),N_geo)),
                      length(theta0),2, byrow=TRUE)

# parameter names
f0 <- function(x) paste0('R0_',x)
f1 <- function(x) paste0('beta_',x)
n_t<- c(sapply(country,f0), sapply(country,f1))

# sd dev for proposal
sigma <- rep(1e-1,length(theta0))

# useful functions
sapply(paste0('Rscript/MCMC/',(list.files('Rscript/MCMC/'))),FUN = source)
##         Rscript/MCMC/adapt.R Rscript/MCMC/adapt_tuning.R
## value   ?                    ?                          
## visible FALSE                FALSE                      
##         Rscript/MCMC/Like1.R Rscript/MCMC/MCMC_full.R
## value   ?                    ?                       
## visible FALSE                FALSE                   
##         Rscript/MCMC/MCMC_iter.R
## value   ?                       
## visible FALSE

run MCMC

#check
rep <- 2e4
# res <- MCMC_iter(iter = rep, theta0 = theta0, s = sigma)
res <- MCMC_full(iter = rep, theta0 = theta0, s = sigma, repli_adapt = 10, within_iter = rep/10)

# save.image(file = 'Rdata/check_mobility.RData')

check convergence

# load(file = 'Rdata/check_mobility.RData')
Acc <- colSums(diff(res$theta)!=0)/(rep-1)
plot(res$logL[,1])

layout(matrix(1:4,2,2))
for (i in 1:length(theta0)){
  plot(res$theta[,i],
       main = paste0(n_t[i],' - ',round(Acc[i]*100)))#,       ylim = prior_theta[i,])
}

Rt daily and Rt_D2

observed mobility

Rt_daily <- array(NA,dim = c(N_days,N_geo,rep))
Rt_D2 <- Rt_daily
m_eff <- Rt_daily
  
for (j in 1:rep){
  R_daily <- Rt_fun(theta = res$theta[j,], x = M )
  Rt_daily[,,j] <- R_daily
  Rt_D2[,,j] <- H %*% R_daily
  
  B <- rep(1,nrow(M)) %*% t(res$theta[j,(2-1)*N_geo+ (1:N_geo)])
  m_eff[,,j] <- 1 + 1/B * log( H %*% exp( - B * (1-M) ) )
}

assumed mobility

n_d <- 1e2
Rt_assumed_mob <- array(NA,dim = c(n_d,N_geo,rep))
x <-  matrix(seq(0,1,length.out = n_d),n_d,N_geo,byrow = FALSE)

for (j in 1:rep){
  R_daily <- Rt_fun(theta = res$theta[j,], x = 1-x )
  Rt_assumed_mob[,,j] <- R_daily
}

summaries

Summary for Rt daily and Rt_D2

# format outputs
temp <- D
temp[,-1] <- NA

results_full_Rt_daily <- list(median_R = temp,
                              low_R = temp,
                              up_R = temp,
                              M = temp)

results_full_Rt_D2 <- list(median_R = temp,
                              low_R = temp,
                              up_R = temp,
                              M_D = temp)

results_meff <- list(median_meff = temp,
                     low_meff = temp,
                     up_meff = temp)

for (i in 1:N_geo){
  
  temp <- apply(Rt_daily[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
  results_full_Rt_daily$median_R[,i+1] <- temp[1,]
  results_full_Rt_daily$low_R[,i+1] <- temp[2,]
  results_full_Rt_daily$up_R[,i+1] <- temp[3,]
  results_full_Rt_daily$M[,i+1] <- M[,i]
  
  temp <- apply(Rt_D2[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
  results_full_Rt_D2$median_R[,i+1] <- temp[1,]
  results_full_Rt_D2$low_R[,i+1] <- temp[2,]
  results_full_Rt_D2$up_R[,i+1] <- temp[3,]
  # results_full_Rt_D2$M_D[,i+1] <- M_D[,i]
  
  temp <- apply(m_eff[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
  results_meff$median_meff[,i+1] <- temp[1,]
  results_meff$low_meff[,i+1] <- temp[2,]
  results_meff$up_meff[,i+1] <- temp[3,]
  
  
}

Summary for Rt assumed mobility

# format outputs
temp <- data.frame(matrix(NA,nrow = n_d, ncol = 1+N_geo))
temp[,1] <- x[,1]
names(temp) <- c('mobility', country)

results_full_Rt_assumed_mob <- list(median_R = temp,
                              low_R = temp,
                              up_R = temp)


for (i in 1:N_geo){
  
  temp <- apply(Rt_assumed_mob[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
  results_full_Rt_assumed_mob$median_R[,i+1] <- temp[1,]
  results_full_Rt_assumed_mob$low_R[,i+1] <- temp[2,]
  results_full_Rt_assumed_mob$up_R[,i+1] <- temp[3,]
  
}

forecast and predictions

short-term

n_pred <- 7
rep_sim <- 1e2
# choose future mobility
M_pred <- rbind(M,matrix(M[nrow(M),],n_pred,N_geo,byrow = TRUE))

# useful functions
sapply(paste0('Rscript/projections/',(list.files('Rscript/projections/'))),FUN = source)
##         Rscript/projections/mobility_prediction.R
## value   ?                                        
## visible FALSE

res_projection <- mobility_prediction(n_pred = n_pred, rep_sim = rep_sim)

temp <- data.frame(dates = seq(1,n_pred) + D$dates[nrow(D)],
                   matrix(NA,n_pred,N_geo))
names(temp) <- c('dates',country)
summary_proj <- list(median = temp,
                           low = temp,
                           up = temp)
for (i in 1:N_geo){
  temp <- apply(res_projection$D_pred[,i,],1,quantile,c(.5,.025,.975))
  summary_proj$median[,i+1] <- temp[1,]
  summary_proj$low[,i+1] <- temp[2,]
  summary_proj$up[,i+1] <- temp[3,]
}

longer-term

first scenario

n_pred <- 60
rep_sim <- 1e2
# choose future mobility
M_pred <- rbind(M,matrix(0.05,n_pred,N_geo,byrow = TRUE))

res_projection_LT1 <- mobility_prediction(n_pred = n_pred, rep_sim = rep_sim)


temp <- data.frame(dates = seq(1,n_pred) + D$dates[nrow(D)],
                   matrix(NA,n_pred,N_geo))
names(temp) <- c('dates',country)
summary_proj_LT1 <- list(median = temp,
                           low = temp,
                           up = temp)
for (i in 1:N_geo){
  temp <- apply(res_projection_LT1$D_pred[,i,],1,quantile,c(.5,.025,.975))
  summary_proj_LT1$median[,i+1] <- temp[1,]
  summary_proj_LT1$low[,i+1] <- temp[2,]
  summary_proj_LT1$up[,i+1] <- temp[3,]
}

Second scenario

n_pred <- 60
rep_sim <- 1e2
# choose future mobility
new_M <- matrix(NA,n_pred,N_geo,byrow = TRUE)
for (i in 1:N_geo){
  x <- 0:(n_pred-1)
  new_M[,i] <- (1-M[nrow(M),i])/(n_pred-1) * x + M[nrow(M),i]
}
M_pred <- rbind(M,new_M)

res_projection_LT2 <- mobility_prediction(n_pred = n_pred, rep_sim = rep_sim)

temp <- data.frame(dates = seq(1,n_pred) + D$dates[nrow(D)],
                   matrix(NA,n_pred,N_geo))
names(temp) <- c('dates',country)
summary_proj_LT2 <- list(median = temp,
                           low = temp,
                           up = temp)
for (i in 1:N_geo){
  temp <- apply(res_projection_LT2$D_pred[,i,],1,quantile,c(.5,.025,.975))
  summary_proj_LT2$median[,i+1] <- temp[1,]
  summary_proj_LT2$low[,i+1] <- temp[2,]
  summary_proj_LT2$up[,i+1] <- temp[3,]
}

Plots

Plotting outputs for each country

for (i in 1:N_geo){
  niceplot(i)
}

plot forecasts

forecast and predictions

short-term

temp <- summary_proj

layout(matrix(1:4,2,2))

for (i in 1:N_geo){
  
 plot_proj(i,threshold = FALSE,log = FALSE)
}

longer-term

first scenario

temp <- summary_proj_LT1

layout(matrix(1:4,2,2))

for (i in 1:N_geo){
  
 plot_proj(i,threshold = TRUE,log = FALSE)
}

Second scenario

temp <- summary_proj_LT2

layout(matrix(1:4,2,2))

for (i in 1:N_geo){
  
 plot_proj(i,threshold = FALSE,log = TRUE)
}

UK

f <- which(country %in% 'United_Kingdom')

layout(matrix(1:4,2,2))

temp <- summary_proj
plot_proj(f,threshold = FALSE,log = FALSE)
temp <- summary_proj_LT1
plot_proj(f,threshold = TRUE,log = FALSE)
temp <- summary_proj_LT2
plot_proj(f,threshold = TRUE,log = TRUE)

containments

when contained??

contained <- data.frame(country = country,
                        median = rep(NA,N_geo),
                        low = rep(NA,N_geo),
                        up = rep(NA,N_geo))

for (i in 1:N_geo){
  
  f <- c(tail( which(results_full_Rt_assumed_mob$median_R[i+1]>1) , 1 ),
    tail(  which(results_full_Rt_assumed_mob$median_R[i+1]>1),1) +1 )
  
  contained$median[i] <- mean(results_full_Rt_assumed_mob$median_R[f,1])
  
    f <- c(tail( which(results_full_Rt_assumed_mob$low_R[i+1]>1) , 1 ),
    tail(  which(results_full_Rt_assumed_mob$low_R[i+1]>1),1) +1 )

  contained$low[i] <- mean(results_full_Rt_assumed_mob$low_R[f,1])
  
    f <- c(tail( which(results_full_Rt_assumed_mob$up_R[i+1]>1) , 1 ),
    tail(  which(results_full_Rt_assumed_mob$up_R[i+1]>1),1) +1 )

  contained$up[i] <- mean(results_full_Rt_assumed_mob$up_R[f,1])
}

write.csv(x = contained,file = 'Rdata/results_contained.csv')

save.image(file = 'Rdata/check_mobility.RData')